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Abstract 

Mathematical  models  of  blood  flow  are  inevitably  embedded  in  models  of  human  ther- 
moregulation because  they  take  the  role  of  the  most  significant  heat  distributor  in  models 
of  the  human  thermal  system  [14,  6].  Models  of  human  thermoregulation  have  a wide 
range  of  applications,  e.g.  for  the  prediction  of  the  impact  of  accidents,  diseases  and  clin- 
ical treatments  (see  [14]  and  the  references  therein).  The  application  of  our  interest  is  the 
prediction  of  the  influence  of  cooling  on  the  heat  distribution  in  premature  infants,  see 
Section  2.  In  Section  3 we  discuss  the  requirements  of  a reliable  thermoregulation  model 
while  the  governing  equation  is  described  in  paragraph  four.  The  employed  blood  flow 
model  is  discussed  within  Section  5.  Section  6 deals  with  numerical  results,  followed  by 
concluding  remarks  in  the  last  paragraph. 


1 Motivation 

Lack  of  oxygen  of  the  fetus  or  newborn  is  known  to  be  an  important  cause  for  injuries  of 
the  developing  brain  [9].  Experimental  studies  have  shown  that  the  neuronal  loss  evolves 
over  several  days  after  such  an  incident  [8].  An  important  factor  influencing  the  degree 
and  distribution  of  neuronal  loss  is  the  cerebral  temperature,  i.e.  lowering  the  cerebral 
temperature  can  prevent  much  damage  [5]. 

The  question  arises,  if  it  is  possible  to  lower  the  cerebral  temperature  of  an  infant  by 

2 - 3 K by  the  manipulation  of  the  environment  inside  an  incubator  while  the  rest  of 
the  body  maintains  a pleasant  temperature.  The  objective  of  this  paper  is  to  discuss  the 
mathematical  measurements  which  can  be  used  to  predict  an  answer  to  that  question 
by  the  use  of  numerical  simulations. 

2 Modeling  the  thermoregulation  of  premature  infants 

The  term  thermoregulation  stands  for  the  measurements  of  the  body  to  hold  a pleasant 
temperature  [4].  Models  for  thermoregulation  consist  of  two  parts:  the  active  and  the 
passive  system  [6].  The  active  system  consists  of  the  regulatory  mechanisms  shivering 
(heat  production  within  the  muscles  attached  to  the  skeleton),  vasomotion  (control  over 
the  degree  of  blood  flow  within  the  skin)  and  sweating  (control  over  the  degree  of  effect- 
iveness of  heat  transfer  between  the  infant  and  the  surrounding  air).  The  passive  system 
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is  the  combination  of  the  physical  human  body  and  the  heat  transfer  in  it  and  at  its 
surface.  The  idea  behind  this  distinction  is  that  the  active  system  has  a controlling  in- 
fluence over  the  passive  system.  Naturally,  only  results  obtained  by  the  complete  model 
can  be  compared  with  available  real  life  data. 

Concerning  premature  infants,  it  is  known  that  shivering  and  sweating  are  not  of 
importance  for  the  modelling  process  [4,  13],  while  vasomotion  should  not  be  of  great 
concern  for  our  special  application  [13].  The  modeling  of  the  passive  system  demands  the 
discretization  of  the  body  and  the  modeling  of  metabolic  heat  production  and  blood  flow. 
We  do  not  consider  phenomena  which  are  related  to  environmental  conditions,  namely 
the  response  to  air  convection,  the  probability  to  gain  or  loose  heat  due  to  radiation  and 
heat  loss  due  to  evaporation  in  dependence  on  pressure,  temperature  and  humidity  of 
the  surrounding  air,  assuming  that  these  are  controllable  by  the  use  of  an  incubator  [13]. 

In  order  to  give  an  answer  to  the  defined  question  by  use  of  numerical  simulations, 
a model  needs  to  deliver  detailed  temperature  profiles  within  the  head  and  a detailed 
resolution  of  the  heat  transfer  processes  in  the  body.  It  should  be  applicable  to  different 
size  neonates  whereby  aspects  like  the  anatomy  and  the  thermal  maturity  have  to  be 
considered.  With  the  exception  of  the  blood  flow  model,  these  aspects  can  be  defined  via 
a suitable  geometry  and  the  use  of  real  life  data  for  spatially  dependent  rates  of  metabolic 
heat  production  within  a numerical  method  [7,  2].  This  also  incorporates  that  existing 
numerical  methods  made  for  the  simulation  of  thermoregulation  of  adults  are  of  no  use 
in  the  given  context  since  studies  have  shown  [3]  that  a detailed  modeling  of  geometry 
and  tissue  composition  is  necessary  in  order  to  obtain  relevant  temperature  profiles.  As 
it  can  be  shown  experimentally  [7,  2]  in  agreement  to  theoretical  discussions  concerning 
thermoregulation  models  of  adults  [6,  14],  the  use  of  a blood  flow  model  greatly  affects 
the  computed  numerical  solutions. 

3 Analysis  of  the  blood  flow  model 

The  bio-heat  equation  derived  by  Pennes  [10]  forms  the  basis  of  the  majority  of  models 
for  human  thermoregulation  in  use  today  [14,  6] . It  describes  the  dissipation  of  heat  in 
a homogeneous,  infinite  tissue  volume.  For  two  spatial  dimensions,  it  can  be  written  in 
the  form 

c(x)p(x)dtT(x,  t)  = div[A(x)VT(x,t)]  + f(x,t).  (3-1) 

Thereby,  the  temperature  T depends  on  the  spatial  variable  x = (xi,  X2)r  as  well 
as  on  time  t.  Furthermore,  A(x),  c(x)  and  p(x)  denote  the  heat  conductivity,  specific 
heat  capacity  and  density  of  the  tissue,  respectively.  The  term  /(x)  can  be  decomposed 
via  /(x,  t)  — Q m (x)  + Qb(x,  t)  into  parts  corresponding  to  metabolic  heat  production 
Qm(x)  and  blood  flow  Qs(x,f). 

As  already  indicated,  the  term  Qm(x)  can  be  defined  by  the  use  of  real  life  data 
[7].  The  formulation  of  the  source  term  due  to  blood  flow  is  based  on  variations  of  the 
following  procedure  [6,  14].  The  idea  is  that  the  body  is  supplied  from  a central  pool  of 
blood  by  the  major  arteries.  Before  the  tissue  is  perfused,  the  temperature  of  the  arterial 
blood  mixes  with  the  temperature  of  venous  blood  flowing  in  adjacent  veins.  After  that, 
the  arterial  blood  exchanges  heat  with  the  tissue  in  the  capillaries  and  becomes  venous 
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blood.  The  venous  blood  is  collected  in  the  major  veins  and  its  temperature  mixes  with 
the  temperature  of  arterial  blood  in  the  adjacent  arteries  before  it  flows  back  into  the 
blood  pool. 

Since  equation  (3.1)  deals  with  the  change  of  thermal  energy  per  unit  volume,  the 
term  QB(x)  takes  the  form 

Qb(x,  t ) - cbPbCCX{x)BF{x)  [TB(t)  - T(x,  t)] , (3.2) 

whereby  TB  (t)  denotes  the  time-dependent  mean  value  of  the  temperature  of  the  blood 
within  the  blood  pool,  we  also  assume  that,  the  specific  density  of  the  blood  pB  and  the 
specific  heat  capacity  of  the  blood  cB  are  constant  variables. 

The  described  modeling  results  in  a differential  equation  for  the  temporal  evolution 
of  the  temperature  within  the  blood  pool,  namely  in 


mBcBdtTB(t)=  [ pBcBCCX{x)BF{x)dx[Tv{t)-TB{t)}.  (3.3) 

Jd 

Thereby,  the  total  blood  mass  mB,  the  time  dependent  mean  value  of  the  temperature 
of  the  venous  blood  Ty{t),  and  locally  defined  tissue-dependent  measures  for  the  blood 
perfusion  BF(x ) and  the  counter-current,  heat  exchange  CCX (x)  are  introduced. 

Equation  (3.3)  shows  that  the  temporal  change  of  the  blood  pool  temperature  is 
proportional  to  the  difference  to  the  temperature  of  the  venous  blood.  The  outlined  idea 
leads  to  the  modeling  of  the  temperature  of  the  venous  blood  as 


fDCCX(x)BF(x)T(x,t)dx 
V[>  fD  CCX (x)BF(x)  d,x 


(3.4) 


which  is  also  usable  when  only  steady  states  are  considered  [7].  The  crucial  terms  in 
the  order  of  importance  are  the  blood  perfusion  BF(x)  and  the  counter  current  heat 
exchange  CCX(x). 

There  is  much  debate  about  the  choice  of  these  functions  in  literature  [14,  6].  This 
debate  arises  because  the  representation  of  blood  circulation  is  substituted  by  a rather 
simple  model  formulation.  The  cure  to  this  disadvantage  is  generally  sought  by  exploring 
more  and  more  detailed  models  of  microstructure,  organs,  etc.,  or  it  is  sought  by  a better 
modeling  of  control  mechanisms  of  the  actice  system  in  the  case  of  adults  [14,  6]. 

The  main  drawback  of  the  described  blood  flow  model  is  given  by  the  blood  pool  idea 
itself.  This  is  up  to  now  to  our  knowledge  not  outlined  in  any  mathematical  description 
of  this  model  within  the  literature  and  can  be  illustrated  as  follows.  Let  a detailed  geo- 
metry be  given  with  a stationary  temperature  distribution  together  with  a homogeneous 
neutral  temperature  at  the  whole  boundary  as  initial  state.  Let  us  assume  that  we  start 
a numerical  computation  where  a selective  cooling  at  the  neck  is  employed.  By  heat  con- 
duction of  the  tissue,  the  effect  of  cooling  computed  with  the  help  of  the  discretization 
of  heat  gradient  and  heat  conductivity  of  the  local  tissue  propagates  into  the  inner  part 
of  the  domain.  Concerning  the  blood  flow,  the  averaging  step  within  (3.4)  captures  the 
local  cooling  effect  which  results  in  a slightly  cooler  average  temperature  of  the  venous 
blood  within  the  whole  domain  than  in  the  initial  state.  Employing  this  value  in  (3.3) 
results  in  a slight  negative  change  of  the  blood  pool  temperature.  Taking  account  of  the 
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evaluation  of  the  source  term  (3.2)  for  the  control  volumes  located  in  the  vicinity  of  the 
neck,  we  notice  that  a strong  cooling  is  locally  equalized  by  the  combination  of  a)  the 
source  term  due  to  blood  flow  which  is  mostly  influenced  by  the  neutral  blood  temper- 
ature in  the  rest  of  the  body  and  b)  of  the  source  term  due  to  metabolic  heat  production 
which  was  not  influenced  at  all  by  the  change  in  the  boundary  temperature.  The  result 
is  that  the  effect  of  a local  cooling  mechanism  is  instantly  distributed  over  the  whole 
domain  while  a weighted  mean  Value  of  the  temperature  over  the  domain  equalizes  local 
cooling  mechanisms.  The  validity  of  this  reasoning  is  verified  by  numerical  results  [7,  2] 
and  by  an  exemplary  result  shown  in  Section  6. 


The  non-local  nature  of  the  described  blood  flow  model  can  directly  be  seen  by 
applying  an  implicit  time  stepping  strategy.  Due  to  the  integration  over  the  whole  com- 
putational domain  in  (3.4),  one  ends  up  with  a fully  occupied  matrix  after  the  usual 
linearization  step  which  was  already  recognized  in  [7]  in  the  context  of  steady  state 
calculations. 


We  now  illuminate  a further  property  of  the  bloodflow  model.  Therefore,  let  the 
abbreviations  a = pbCb,  fl  = fD  Kb(x)B(x)  dx  and  7 = pB/mB  hold.  A straightforward 
computation  gives 


(3.5) 


Note  that  a,  (3  and  7 are  positive  constants.  Consider  a steady  state  situation  as  initial 
state,  i.e.  Tb  = Ty  holds.  If  the  body  is  heated,  the  temperature  within  the  body 
increases  and  so  Ty  will  increase.  This  has  the  effect  that  the  bloodpool  temperature  Tb 
will  increase  in  the  near  future,  i.e.  t'b  (t)  > 0.  We  now  investigate  the  net  effect  of  the 
bloodflow.  Integration  of  the  source  over  the  computational  domain  D results  in 


b 


QB{x,t)  dx  = a 


flTB(t)  — [ Kb  (x)B(x)T(x,  t)  dx 
Jd 


(3.5)  ad 

= “ dtTB{t)- 


When  employing  T'B{t)  > 0 we  see  that  the  total  of  all  sources  in  the  body  is  negative, 
i.e.  while  the  blood  in  the  bloodpool  cools  the  increasingly  warm  body  in  the  mean  if 
the  body  is  exposed  to  heat,  it  also  takes  over  heat  from  it.  The  bloodpool  and  the  body 
are  to  be  seen  as  two  separate  systems  which  are  connected  via  heat  fluxes  and  so  one 
can  consider  the  bloodpool  as  a regulator. 


4 Numerical  method  and  experiments 

The  following  numerical  approximation  of  the  unsteady  bio-heat  equation  (3.1)  repres- 
ents a convenient  extension  of  the  finite  volume  method  developed  in  [7],  which  has  been 
proven  to  be  a robust,  accurate  and  reliable  algorithm  in  the  context  of  steady  state 
temperature  distributions.  However,  finite  volume  schemes  are  categorically  based  on 
the  integral  form  of  the  governing  equation.  In  order  to  apply  Gauss’s  integral  theorem 
it  is  neccessary  to  write  the  equation  in  divergence  form.  Therefore,  we  introduce  the 
auxiliary  variable  k(x)  = p(x)c(x)  and  the  auxiliary  temperature  T(x,t)  = k(x)T(x,t ) 
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into  the  governing  equation  and  consequently  the  bio-heat  equation  (3.1)  writes 


J T(x,t)dx  = ^jvr(x,t)  - vfc(x)  • n(x)  ds  + J f(x,  t) 


d_ 
dt 


I dx 
(4.1) 

for  all  control  volumes  a C D,  see  [2],  In  order  to  solve  equation  (4.1)  numerically, 
the  space  part  D is  decomposed  into  a finite  number  of  sub-domains.  We  start  from  an 


Fig.  1.  General  form  of  a control  volume  of  the  triangulation  (left)  and  its  boundary 
(right). 


arbitrary  conforming  triangulation  T>h  of  the  domain  D which  is  called  the  primary  mesh 
and  consisting  of  finitely  many  triangles  T>i  and  the  corresponding  nodes  are  abbreviated 
by  x,  € D.  Based  on  the  triangulation  a discrete  control  volume  ct,  is  defined  as  the  open 
set  of  R2  including  the  node  x,  and  bounded  by  straight  lines  which  are  determined  by 
the  connection  of  the  midpoints  of  the  edges  of  the  corresponding  triangles  Vj  (i.e. 
Xj  € dVj)  and  their  barycentre  (see  Figure  1).  The  union  Bh  of  all  boxes  is  called  the 
secondary  mesh.  A finite  volume  method  represents  a discretization  of  the  evolutionary 
equation  (4.1)  for  cell  averages  defined  by  ( MT ) (f)|ff  = (1/M)  faT(x,t)  dx,  where  \a\ 
denotes  the  volume  of  the  box  a.  With  respect  to  the  secondary  mesh  Bh  we  can  write 
the  integral  form  (4.1)  as 


d_ 

dt 


(Afr)(,)l'.  = ra{/sJ®vr(x-,)-^wiav‘H'n(x)ds 


+ 


f QB(x,t)dx+  f QM(x)dx\ , V<ji  e Bh. 

J <7i  J <7i  ) 


(4.2) 


Corresponding  to  a finite  element  method  the  evaluation  of  the  boundary  integral  is  per- 
formed by  using  a piecewise  const  ant  distribution  of  the  heat  coefficient,  A and  a piecewise 
linear  distribution  of  the  auxilary  temperature  T.  with  respect  to  the  triangles  of  the 
trianglutaion  used.  Note  that  the  source  term  remains  unchanged  and  the  calculation  is 
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given  by 


and 


dx  = \ai\QM(x.i) 


W^csPBCCX^BFixi)  [TB(t)  - T(xj,t)] . 


The  computation  of  the  blood  pool  temperature  is  directly  performed  by  an  explicit  time 
discretization  of  e_quation  (3.3).  Thereby,  the  temperature  of  the  venous  blood  is  given 
by  equation  (3.4). 

It  is  remarkable  that  the  method  degenerates  to  the  scheme  presented  in  [7]  in  the 
context  of  a steady  state  solution  and  therefore  the  excellent  properties  like  the  dis- 
crete min-max  principle  are  maintained  in  such  a situation.  Due  to  the  space  available 


Fig.  2.  Primary  mesh  and  tissue  layers  in  the  head  region. 

we  restrict  ourself  to  the  consideration  of  steady  state  calculations  using  the  described 
method.  Thereby,  we  distinguish  layers  of  skin,  fat,  bone  and  kernel  by  different  rates 
of  metabolism,  specific  heat  capacity  and  blood  perfusion  associated  with  the  regions 
depicted  in  Figure  2.  As  boundary  conditions  we  employ  a comfortable  boundary  tem- 
perature of  309.15  K at  head,  back,  legs,  and  belly  while  we  set  299.15  A"  at  the  neck, 
i.e.  we  selectively  cool  the  neck.  In  reality,  this  corresponds  to  the  situation  where  the 
infant  is  wearing  a water-filled  collar  with  the  purpose  of  cooling  the  blood  flowing  into 
the  brain  through  the  arteries  adjacent  to  the  skin. 

In  Figure  3 (a)  we  can  see  the  temperature  distribution  in  the  two-dimensional  dis- 
cretized idealization  of  the  body  of  a premature  infant.  Thereby,  no  blood  flow  and  no 
metabolic  heat  production  is  applied,  so  that  the  depicted  distribution  of  heat  is  only 
influenced  by  the  heat  conductivity  of  the  employed  tissues.  The  situation  where  tissue 
dependent  metabolic  heat  production  is  taken  into  account  is  shown  in  Figure  3 (b). 
Note  that  the  heat  sources  visualized  within  the  picture  not  only  have  local  effects,  they 
also  influence  the  mean  value  of  the  temperature  of  the  blood  pool.  Within  Figure  3 (c), 
blood  flow  is  additionally  given. 
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It  is  evident  that  the  blood  flow  has  the  effect  outlined  in  Section  5.  Especially,  the 
numerical  solution  incorporates  no  hint  of  the  fact,  that  in  reality  there  is  a transport 
of  cool  blood  to  the  brain  and  also  a transport  of  blood  by  the  veins  coming  from  the 
brain. 


299.15  301.1  303.1  305.1  307.0  309.0  310.3 

Temperature  in  [K] 

Fig.  3.  Comparison  of  steady  state  situations  (a)  only  with  heat  conduction  (b)  with 
heat  conduction  and  metabolic  heat  production  and  (c)  with  blood  flow  additionally 
taken  into  account  (from  top  to  bottom). 


5 Concluding  remarks 

The  range  of  applicability  of  the  described  blood  flow  model  is  restricted  to  situations 
where  it  makes  sense  to  employ  a mean  value  of  the  whole  blood,  e.g.  if  the  whole  body  is 
exposed  for  a longer  time  to  the  same  temperature.  For  a clinical  application  where  the 
effects  of  local  cooling  or  heating  have  to  be  studied,  caution  is  required  when  dealing 
with  the  results  achieved  by  employing  variations  of  the  described  model. 
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